function [S, info] = slgraphupd(G, S0, cf0, varargin)
%SLGRAPGUPD Performs Graph-based state updating through propagation
%
% [ Syntax ]
%   - S = slgraphupd(G, S0, cf0, ...)
%   - [S, info] = slgraphupd(G, S0, cf0, ...)
%
% [ Arguments ]
%   - G:        The graph adjacency matrix (n x n)
%   - S0:       The initial states (n x d)
%   - cf0:      The confidences on initial states (1 x 1 or 1 x n)
%   - info:     The struct of the procedural information
%               It has the following fields:
%               \{:
%                   - method:   the method used to update the states
%                   - niters:   the number of iterations
%                   - converged: whether the solution converges
%               \:}
%
% [ Description ]
%   - S = slgraphupd(G, S0, cf0, ...) performs graph-based updating through
%     local propagation. 
%
%     Suppose there are n nodes in the graph, and each node is associated
%     with a d-dimensional state vector. Then G should be an n x n matrix,
%     S0 should be an n x d matrix with each row giving the states of a
%     vector, and cf0 is a 1 x n row vector giving the confidence values on
%     initial states. If all source vectors have the same initial
%     confidences, cf0 can also be a scalar.
%
%     The updating can be considered as a procedure that repeatedly applies
%     the following rule
%     $ S(i,:) <-- (1 - cf0(i)) * \sum_j G0(i,j) * S(j,:) + cf0(i) * S0(i,:) $
%
%     It can be seen that at each step, the new state is set to a weighted 
%     combination of the initial state and the neighboring states.
%
%     In implementation, there are two ways to accomplish the updating. One
%     is to directly compute the close-form converged solution which
%     involves the spectral analysis of the graph; the other way is to
%     follow the iterative updating formula. 
%
%     Typically, direct computation of the close-form solution is more
%     accurate and convenient for small graphs. However, it becomes
%     instable (due to singularity of the adjancency matrix) and
%     inefficient when the problem scale gets larger. For large problem, it
%     is preferrable to follow the iterative procedure, and tune the
%     trade-off between accuracy and run-time by adjusting the convergence
%     criteria.
%
%     You can specify the following options for the updating implementation
%     \{:
%         - method:         The method to implement the updating.
%                           (default = 'csol')
%                           \{:
%                               - csol:  directly compute the close-form
%                                        solution
%                               - iter:  iterative updating
%                           \:}
%
%         - maxiter:        The maximum number of iterations. (default = 100)
%                           The option only takes effect for iter method.
%
%         - tol:            The tolerance of change of convergence.
%                           (default = 1e-6)
%                           It is measured by the maximum value difference.
%                           The option only takes effect for iter method.
%
%         - convpred:       The function to determine whether the procedure
%                           converges. If is invoked by the syntax as
%                              $ f(S_old, S_new) $
%                           Here, S_old and S_new are respectively the
%                           state vector matrix before/after the current
%                           iteration.                              
%                           The option only takes effect for iter method.
%
%         - verbose:        Whether to display iteration information.
%                           (default = false)
%                           This option only takes effect for iter method.
%     \:}
%
%   - [S, info] = slgraphupd(G, S0, cf0, ...) additionally returns the
%     relevant information of the updating procedure.
%
% [ Remarks ]
%   - Here briefly explains the mathematical procedure of the constraint on
%     the input graph.
%
%     The aforementioned updating rule can be concisely written in form of
%     matrices as follows
%       $ S <- (I - D_cf) * G0 * S + D_cf * S0 $
%     Here, D_cf is diag(cf0). 
%
%     It can be proved that a converged solution can be derived if and only
%     if the spectral radius of the updating matrix U = (I - D_cf) * G0 is
%     less than 1. 
%
%     A typical construction is to create a graph G with spectral radius
%     equal to or less than 1, and let all confidence values set in the
%     range [0, 1]. For each connected component of G with spectral radius 
%     equaling 1, it should be guaranteed that at least one node in this
%     component has positive initial confidence.
%
%   - If convpred is specified, then convpred will be used to judge the
%     convergence, otherwise convergence will be judged by comparing the
%     maximum value change with tol.
%
% [ History ]
%   - Created by Dahua Lin, on Jul 8, 2007
%

%% parse and verify input arguments

% G, S0, cf0
assert(isnumeric(G) && ndims(G) == 2 && size(G, 1) == size(G, 2), ... 
    'sltoolbox:slgraphupd:invalidarg', 'The graph G should be a numeric square matrix.');
assert(isnumeric(S0) && ndims(S0) == 2, 'sltoolbox:slgraphupd:invalidarg', ...
    'The S0 should be a 2D numeric matrix.');
assert(isnumeric(cf0) && ndims(cf0) == 2 && size(cf0, 1) == 1, ...
    'sltoolbox:slgraphupd:invalidarg', 'The cf0 should be a numeric scalar or row vector.');

n = size(G, 1);
assert(size(S0, 1) == n, 'sltoolbox:slgraphupd:sizmismatch', ...
    'The number of rows in S0 should equal the number of nodes.');
assert(isscalar(cf0) || size(cf0, 2) == n, 'sltoolbox:slgraphupd:sizmismatch', ...
    'The cf0 should be either a scalar or a row vector of length n.');

% options

opts = struct( ...
    'method', 'csol', ...
    'maxiter', 100, ...
    'tol', 1e-6, ...
    'convpred', [], ...
    'verbose', false);
opts = setopts(opts, varargin{:});

assert(ischar(opts.method), 'sltoolbox:slgraphupd:invalidarg', ...
    'The option method should be a string.');
use_iter = find(strcmp(opts.method, {'csol', 'iter'})) == 2;
assert(~isempty(use_iter), 'sltoolbox:slgraphupd:invalidarg', ...
    'The option method should be eieht ''csol'' or ''iter''.');


%% main delegate

if ~use_iter
    if nargout < 2
        S = gupd_csol(G, S0, cf0);
    else
        [S, info] = gupd_csol(G, S0, cf0);
    end
else
    if nargout < 2
        S = gupd_iter(G, S0, cf0, opts);
    else
        [S, info] = gupd_csol(G, S0, cf0);
    end
end


%% The close-form solution method

function [S, info] = gupd_csol(G, S0, cf0)

[U, B] = calc_ub(G, S0, cf0);

A = eye(size(U)) - U;
S = A \ B;

if nargout >= 2
    info = struct( ...
        'method', 'csol', ...
        'niters', 0, ...
        'converged', true);
end


%% The iterative procedure method

function [S, info] = gupd_iter(G, S0, cf0, opts)

% check arguments
mxit = opts.maxiter;
assert(isnumeric(mxit) && isscalar(mxit) && mxit > 0 && mxit == fix(mxit), ...
    'sltoolbox:slgraphupd:invalidopt', 'The option maxiter should be a positive integer scalar.');

tol = opts.tol;
assert(isnumeric(tol) && isscalar(tol) && tol > 0, ...
    'sltoolbox:slgraphupd:invalidopt', 'The option tol should be a positive numeric scalar.');

if isempty(opts.convpred)
    fconv = @(s0, s1) max(abs(s1(:) - s0(:))) <= tol;
else
    assert(isa(opts.convpred, 'function_handle'), ...
        'sltoolbox:slgraphupd:invalidopt', 'The option convpred should be empty or a function handle.');
    fconv = opts.convpred;
end

% iterate updating
nit = 0;
S = S0;
converged = false;

[U, B] = calc_ub(G, S0, cf0);

if opts.verbose
    fprintf('Start iterative updating.\n');
end

while ~converged && nit < mxit     
        
    nit = nit + 1;
    S_old = S;
    
    S = U * S + B;
    
    converged = fconv(S_old, S);
    
    if opts.verbose
        chx = max(abs(S(:) - S_old(:)));                
        fprintf('Iter %5d: max change = %g\n', nit, chx);
        if converged
            fprintf('The iterative updating converges.\n');
        end
    end
end

if nargout >= 2
    info = struct( ...
        'method', 'iter', ...
        'niters', nit, ...
        'converged', converged);
end



%% Auxiliary computation function 


function [U, B] = calc_ub(G, S0, cf0)

if isscalar(cf0)
    U = (1 - cf0) * G;
    B = cf0 * S0;
else
    U = bsxfun(@times, G, 1 - cf0(:));
    B = bsxfun(@times, S0, cf0(:));
end

